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The critical exponents for a class of one-dimensional models of interface depinning in disordered 
media can be calculated through a mapping onto directed percolation (DP). In higher dimensions 
these models give rise to directed surfaces, which do not belong to the directed percolation uni- 
versality class. We formulate a scaling theory of directed surfaces, and calculate critical exponents 
numerically, using a cellular automaton that locates the directed surfaces without making reference 
to the dynamics of the underlying interface growth models. 
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The last decade's intense theoretical, numerical, and 
experimental interest in the growth and roughening of 
interfaces has been fueled in part by the interdisciplinary 
aspects of the subject [Q. Applications include fluid- fluid 
displacement Q , imbibition in porous media || , and the 
motion of flux lines in superconductors 0. In such sys- 
tems, an interface advances through a disordered (typi- 
cally porous) medium, under the influence of an external 
driving force, F. There exists a critical value, F c , of this 
force, such that for F < F c the interface is pinned by 
the disorder, while for F > F c it moves with a constant 
velocity v. In the vicinity of the depinning transition 
(F = F c ), v ~ f, where / = (F - F c )/F c . The cor- 
relation length (£) characterizing the size of the pinned 
regions in the plane of the interface diverges at F c as 
£ ~ f~ v ■ For F > F c the width of the interface in steady 
state varies with the system size L as w(L) ~ L a . Here 
9, v and a are the velocity, correlation length, and rough- 
ness exponents, respectively. 

In this paper we consider "directed percolation de- 
pinning " (DPD) models |,|, which are believed to de- 
scribe the depinning transitions in a variety of systems, 
among them interfaces described by the Kardar-Parisi- 
Zhang equation 01 with quenched noise and imbibi- 
tion experiments |3] . The DPD model is defined on a 
d-dimensional lattice @], with periodic boundary condi- 
tions in the (d— 1) interface dimensions and open bound- 
ary conditions in the dimension corresponding to the di- 
rection of motion of the interface (the z direction, say). 
Sites are randomly occupied by impurities with proba- 
bility p, and a fluid is imagined to push its way upward 
from below. At each time we randomly choose one of 
the impurity-free "dry" sites neighboring the interface 
that separates the wet and dry regions. The interface 
advances by "wetting" the chosen site, and any site be- 
low it in the same column (along the z direction) ||. In 
this model only continuous, unbroken surfaces of impuri- 
ties that cover the entire (d-l)-dimensional cross section 
of the lattice and contain no overhangs can successfully 
block the interface . In 2d, this blocking surface is a 
collection of Id lines, which can be viewed as the back- 
bone H of the infinite cluster of a (l-fl)-dimensional di- 
rected percolation (DP) problem that starts from one side 
of the lattice and percolates towards the other. Thus the 
scaling exponents (a, 9, v) for the DPD model in 2d can 



be obtained from the mapping to the (l+l)-dimensional 
DP problem In particular, the correlation length of 
DPD can be identified with the longitudinal correlation 
length of DP, giving v = vn DP w 1.73. Similarly, one ob- 



tains the exponents a and 9 as a = v± /v\\ « 0.63, 



and 9 = v\\ UF - i/ x « 0.63, where u x r ~ blO is the 
transverse correlation length exponent of DP. 

While the DP theory correctly predicts all relevant ex- 
ponents for 2d DPD models, the mapping fails in higher 
dimensions. There is a simple topological reason for this: 
In (l+l)-dimensions, the directed percolation backbone 
is a collection of Id lines, capable of blocking the mo- 
tion of the one-dimensional interface in 2d DPD. In 3d, 
however, only an unbroken 2d surface of pinning sites 
without overhangs can block the advance of the 2d inter- 
face. The collection of such surfaces forming the back- 
bone of a blocking cluster in 2d (or in any other dimen- 
sion) is referred to as a directed surface [JlOj . Numerical 
measurements of the critical exponents characterizing the 
depinning transition of DPD models for d > 2 jll]] con- 
firm that directed surfaces belong to a universality class 
different from DP. 

In this paper we take several steps toward a descrip- 
tion of the DPD depinning transition for d > 2 in terms 
of directed surfaces. First we introduce a deterministic 
two-state cellular automaton (CA) that finds the directed 
surface for a system with an arbitrary distribution of im- 
purities, and in arbitrary dimension. The CA produces 
this surface without making reference to the dynamics of 
the underlying DPD interface growth models. Focusing 
on the directed surfaces in this way allows us to formulate 
a scaling theory of the transition, and thus describe this 
problem using the standard formalism of critical phenom- 
ena. Moreover, the CA allows us to measure numerically 
exponents not available from the DPD models. Using the 
derived scaling laws we thereby obtain a complete set of 
scaling exponents characterizing the static properties of 
directed surfaces. 

Cellular automata and directed surfaces — We first de- 
fine the CA rules, and then explain why they generate the 
desired directed surfaces. Consider a square 2d lattice in 
the {x, z) plane. At time t = each site (x, z) is either in- 
dependently occupied by an impurity {sq(x, z) = 1) with 
probability p, or is empty (sq(x, z) = 0). The CA rule 
is defined as follows: s t +i{x, z) = 1 if the following three 
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conditions are simultaneously satisfied: 
s t (x, z) = 1; 

s t (x — 1, z — 1) + St(x — 1, z) + st(x — 1, z + 1) > 0; 

s t (z + 1, z - 1) + s t (a; + 1,4 + + 1, z + 1) > 0. 

If any of these is not satisfied, then s t +i(x, z) = 0. 
Boundary conditions are periodic in x and open in z. The 
CA rule leaves untouched any occupied site that locally 
belongs to an unbroken Id path roughly perpendicular to 
the z direction. If an occupied site (x, z) has an occupied 
neighbor or next near neighbor in both the (x — l)st and 
(x+l)st columns (see Fig. 0(a)), then it is part of such a 
path. However, if the occupied site is at the tip of a dan- 
gling branch (e.g., the site B in Fig |l](a)), then at the next 
time step it will be removed, since it is missing a neigh- 
bor in one of the adjacent columns. Once B is removed, 
however, A is left without a neighbor in one adjacent col- 
umn, and so is itself removed at the following time step. 
In this fashion, all dangling branches and isolated clusters 
of impurity sites are systematically eliminated. Thus un- 
der successive applications of the rule, all impurities that 
do not belong to the backbone of the infinite cluster of 
the equivalent (l+l)-dimensional DP problem (e.g., the 
shaded path in Fig. |l|(a)), wherein z and x respectively 
correspond to space and time, ultimately disappear. If 
such an infinite DP cluster does not exist, then the fixed 
point of the CA has all sites empty. Fig. |l|(b) shows the 
result of applying this rule in a system of size L = 100 
with p = 0.55. 

The generalization to higher dimensions is straightfor- 
ward. We discuss only the 3d case, defined on a cubic 
lattice with axes labeled (x, y, z); the extension to d > 3 
is obvious. In 3d, s t +\{x,y, z) = 1 if the following five 
conditions are simultaneously satisfied: 

s t {x,y,z) = 1; 

s t (x-l,y,z-l)+s t (x-l,y, z) + s t (x-l,y, z+1) > 0; 

s t (x+l,y, z-l)+st(x+l,y, z) + s t (x+l,y, z+1) > 0; 

st(x,y — l,z-l) + s t {x,y-l,z) + s t (x,y — I, z+1) > 0; 

st(x,y+l,z-l) + s t (x,y + l,z) + st(x,y + l,z+l) > 0. 

If at least one of these conditions fails, then 
s t+ i(x, y, z) — 0. Periodic and open boundary conditions 
apply in the {x,y) and z directions, respectively. 

As in 2d, for d > 3 the CA rule removes any site that 
does not belong to a locally continuous (d-l)-dimensional 
surface perpendicular to the z direction. It retains only 
unbroken surfaces, without overhangs, capable of block- 
ing the advance of an interface in the associated DPD 
model. It therefore locates exactly all directed surfaces in 
arbitrary dimension. Thus one expects the critical expo- 
nents measured for the DPD model and for the directed 
surfaces generated by the CA to coincide. This conclu- 
sion is supported by the numerical results presented be- 
low. 

Scaling exponents for directed surfaces- The final state 
(fixed point) of the CA depends on p. For small p 
there are no directed surfaces in the system, so the aver- 
age density of the final state, p(p) — J2 r s ( r )' wnere 
r = (ii, ii, --id), is zero. Increasing p, one reaches a criti- 
cal value, p c , so that for p > p c , p{p) is nonzero. In anal- 



ogy with percolation, one expects that p(p) ~ (p — p c Y ■ 
The exponent (3 is not known, and has not previously 
been measured numerically, except in 2d, where the di- 
rected surface is the backbone of the infinite cluster of 
(l+l)-dimensional DP. Since the exponent (3 for this 
backbone is known to be 2(3 DP |§, where f3 DP ~ 0.28 
in (1+1) dimensions, we have /? « 0.56 in 2d. 

To characterize further the directed surfaces for ar- 
bitrary d, it is helpful to define parallel (to the aver- 
age orientation of the surface) and perpendicular cor- 
relation lengths, £|| and respectively. Near the de- 
pinning transition they diverge as £y ~ (p — p c ) and 

U~(P-Pc)- VX - 

Finally, one can define exponents associated with the 
distribution of the sizes of voids or holes in the directed 
surface, i.e., empty regions totally surrounded by sur- 
face sites ]Iiyr3[ |. The probability distribution function, 
P v (v), for the void volumes v, and Pi(s) for linear void 
sizes s in the ith direction, are expected to behave alge- 
braically: P v (v) ~ v~ Tv , and, in 3d, e.g., P x . y (s) ~ s~ r n , 
and P z (s) ~ s~ T± . As recent work by Huber et al. shows 
flTij , in 2d the void exponents r„, ry, and t± can be re- 
lated to the exponents characterizing the size distribu- 
tion of avalanches associated with the dynamics of the 
DPD model Jll]] or its self-organized version, the "self- 
organized depinning" model |l3| , |l5| ] . 

Scaling theory- The first advantage of describing di- 
rected surfaces in static terms similar to those fruitfully 
employed for the percolation problem is that we can use 
the standard scaling arguments familiar from critical phe- 
nomena to characterize the correlations and the fractal 
nature of those surfaces Jl6| . As a first step, we derive 
a set of scaling relations for the critical exponents of di- 
rected surfaces. Imagine a coarse-grained density field 
ip(x) for the surface. Under rescaling of distances in the 
directions parallel and perpendicular to the surface via 
X[| = and x± = b a x' ± , the field tp and the distance 
A = p — p c from the critical point are assumed to rescale 
according to ip(x\\ , x±_) = b x tp'(x^ , x' ± ) and A = b^ 1 ^ A', 
respectively. Here b(> 1) is the length rescaling factor, 
and a, and \ are critical exponents. 

The rescaling of lengths implies that the correlation 
length exponents are given by v» = and i>±_ = uQ: De- 
noting the average density of the surface, < ip(x) >, by 
Al, it follows that M vanishes like A' 3 as A — > + , with 
P = ~XV\\- 

To express x m terms of the exponents for other cor- 
relation functions, note that the scaling relations above 
imply 

G(x, A) = A- 2 *"" h(x\\ A"n , x±A v± ) (1) 

where G(x, A) =< ^{x)ip{Q) >. For A > 0, G has two 
parts, a disconnected piece equal to M 2 , and a connected 
piece G c (x, A) = G(x, A) — M 2 . G(x, A) is proportional 
to the steady-state probability that both sites and x 
belong to the directed surface, i.e., to the product of the 
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probability of site belonging to the directed surface 
and the conditional probability that x belongs to the di- 
rected surface, given that does. This makes it clear that 
G(x, A) vanishes like M as A — > 0, whereupon, close to 

p c , G c (x\\ t ±) ~ Mi. ^ I/ "'" L . Using the standard notation 
of critical phenomena, we define the exponents rj\\ and 



i?_L via G c {x\\. 
laws 



.) - Mi 
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leading to the scaling 



/3 = 2/||,x(rf-2 + r?||,x) 



(2) 



The exponents 7711 and r]± are straightforwardly related 
to the fractal dimensions |l8| ], and of the directed 
surface in the || and _L directions by [|l9| D\\ = 1 — T)\\ and 
D± = 3 — d — rj±. Moreover, the overall fractal dimension, 
D, defined by the total number of surface points within 
a distance R of a given point on the surface growing like 
R D , is given by D = 2 — rj± = d — (3/v±- 

Finally, for voids, we have the scaling relations rii = 
l + (T v -l)(d-l + a), andri = 1 + (r - l)(d- l + a)/a. 
In 2d the void exponents can be related to V\\ , and /3, 
through the formula 1 14 1 t v — \ = (v\\+v_i_ — 2/3)/(v\i+v±). 

Numerical results- In 2d all the exponents defined here 
are available from the DP analogy. For higher d, vn, i/±, 
and hence a have been obtained from DPD simulations 
|pr| . Others, such as (3, rjux and Du± are difficult to 
compute from DPD, and so have not been measured. 
The CA representation makes numerical determination 
of these rather straightforward. Moreover, the scaling 
law (Q) shows that, given 1/11 and one additional in- 
dependent exponent suffices to fix all the others, except 
the void exponents, which for d > 2 have so far not been 
related to the others. 

In principle, the most straightforward exponent to 
measure numerically is (3, obtained from plotting the den- 
sity of the final state of the CA as a function of p. In 2d 
the data show a fairly unambiguous scaling regime, yield- 
ing the value (3 « 0.5, consistent with the known value 
f3 BB = 2f3 DP r* 0.56. In 3d, near p c , p shows a rather 
precipitous jump that sharpens with increasing sample 
size, suggestive of either a first order phase transition or 
a very small exponent (3. The data for p(jp) are insuffi- 
cient for one to choose between these possibilities. We 
therefore resort to determining (3 by measuring D|| and 
D±, and using the scaling relations D» = d — 1 — (3/v\\ 
and D±_ = 1 — /3/i/±, and the values = 1.18 ±0.10 and 
v± = 0.57 ±0.05, taken from DPD simulations Q. The 
values for (3 thus obtained, together with our values of 
Du and D±, are listed in Table |. Our numerical values 
for the void exponents are also given in the table. The 
existence of a fairly clear scaling regime for Dm with a 
value less than 2 argues against the occurrence of a first- 
order phase transition (as, of course, does the continuous 
phase transition observed in the 3d DPD model). The 
scaling observed in the data of Fig. 2b (spanning, for P v , 
more than two orders of magnitude in void sizes), also 
strongly supports the inference of a second-order transi- 
tion. The 2d value for t v is consistent with the prediction 



t v = 1.80 by Huber et al. [[Uj. Using the 3d numerical 
values, we estimate from the scaling relations above the 
value a ~ 0.48, in good agreement with the value ob- 
tained from the DPD model. This provides further con- 
firmation that the directed surface generated by our CA 
indeed belongs in the same universality class as DPD. 

One of the major benefits of the CA approach intro- 
duced in this paper is that it replaces the dynamic models 
used to study the properties of directed surfaces with a 
static picture. Dynamic models cannot capture all di- 
rected surfaces existing in the system. Models with ran- 
dom updating, such as DPD, remove small patches of the 
directed surface in an uncontrolled fashion during large 
avalanches. The SOD models Jll| also systematically 
eliminate branches of directed surfaces. By contrast, the 
CA is the first algorithm that identifies all underlying di- 
rected surfaces, which in turn determine static exponents 
such as v\\ , vj_, and a, some of which (such as (3), are not 
accessible from dynamic models. The CA also allows di- 
rect calculation of the void-size exponents, from which, 
at least in 2d, the avalanche exponents of the dynamic 
models can be derived ||l4j| . 

We thank Sid Redncr for helpful discussions. 
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FIG. 1. (a) Schematic representation of the action of the 
CA rule. The final state of the CA contains only the back- 
bone of the underlying DP cluster (solid circles). Open circles 
represent dangling branches or isolated clusters that will even- 
tually be eliminated, (b) Directed surfaces produced by the 
CA rule in 2d with p = 0.55 and L = 100. 

FIG. 2. Void size distributions determined numerically us- 
ing the CA rule, (a) 2d: The upper, middle and lower curves 
correspond to P v (v), P x {s) and P z (s), respectively. We used 
p = 0.54, L = 1000 x 1000, and averaged over 250 runs, 
(b) 3d: The upper and lower curves correspond to P v (v), 
and Pz(s), respectively. In the middle we have two curves 
superposed, corresponding to P x ,y(s). We used p — 0.74, 
L = 100 x 100 x 100, and averaged over 1500 runs. 



TABLE I. Exponents obtained from numerical simulations 
using the CA rule, and from scaling relations. 



Exponent 


2d 


3d 


P 


0.5 ±0.07 


0.1 ±0.02 




0.71 ± 0.03 


1.9 ±0.05 


D ± 


0.54 ± 0.05 


0.8 ±0.1 


T v 


1.7 ±0.1 


2.2 ±0.1 


Tx 


2.2 ±0.1 


3.9 ±0.2 


Tz 


2.8 ±0.2 


7.0 ±0.4 
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